Séries temporelles

Projet 1 : detection des dates de rupture dans le ROA

Projet 2 : test de racines unitaires et predictions sur le ROA

Mettre une intro : type de la série, info sur le pays sur la période…

Afficher les crises sur un chronogramme. Pas de monthplot ou decompose car données annuelles

Calculer les stats descriptives de la série.

Tests de racines unitaires (les 4)

Dans le corps du texte, on ne met que la bonne spécification (la dernière de l’arbre de décision du cours) du test de racines unitaires, et les autres en annexes. Il faut justifier (juste une phrase) tous les résultats, même ceux en annexes. Mettre les équations des modèles.

Mettre l’équation du modèle estimé

Interpréter une statistique de test, au moins une fois : H0, H1, la formule, et la distribution sous H0.

Rendre la série stationnaire, une fois que c’est le cas : ACF, PACF, LUNGBOX

Pour chaque graphique : commenter le graph forcément

Si autocorrélation : modéliser cette dernière avec EACF. Cette dernière nous donne les valeurs de départ de p et de q et donc on fait ARMA(p,q). A la fin, tous les coefficients du modèle doivent être significatifs. S’il existe un modèle mieux que le notre, on perd 1 point.

Si annexe remplie, point bonus.

Tous les modèles, même annexe, doivent être commentés.

Pour les résidus

Ensuite on vérifie que les aléas sont des bruits blancs. E() = 0. Si on rejette H0, cela veut dire que l’on a enlevé à tort l’intercept et qu’il faut le rajouter)

Test de ljung box pour vérifier qu’il n’y a plus d’autocorrélation dans les résidus. Si on rejette H0, alors on ré estime le modèle ARMA(p,q) avec des p et q plus grands.

Test ARCH sur l’homoscédasticité.

Si les 3 H0 sont vrais, on peut faire de la prévision : d’abord sur la série transformée, puis sur la série originale.

Tester jacques berra pour la normalité, sinon, test agostino et anscombe sur les résidus + histogramme (on cherche à savoir si le rejet de l’hypothèse associée à jarque berra est dû à la skewness ou au kurtosis)

Comparer des modèles

• tous les coefs doivent être significatifs
• les aléas sont des bruits blancs

On choisit celui qui minimise le BIC (possibilité de choisir d’autre critère, mais il faut justifier)

Ne jamais parlé de zlag1 ou difflag1, mais toujours leur nom dans le cours (beta0 etc)

Projet 2 : si les deux processus sont DS, alors la régression est fallacieuse (solution : cointégration et mce)

Pour pouvoir faire les MCO, il faut que toutes les variables utilisées soient issues d’un pgd stationnaire. Dans le cas de PGD ds, la stat t ne nuit pas une student.

Contexte

On cherche à savoir si le PGD est stationnaire. Si oui, alors :

• \(E(y)\) ne dépend pas du temps
  
• \(Var(y)\) finie
  
• Autocovariance entre \(Y_t\) et \(Y_{t-h}\)

S’il n’est pas stationnaire, on change la série pour qu’elle le devienne et qu’on puisse faire notre estimation. Pour savoir s’il est stationnaire, on fait des tests de racines unitaires.

Tendance TS : \(X_t = \alpha_0 + \alpha_1 X_{t-1} + \varepsilon_t\)

Tendance DS : \(X_t = X_{t-1} + \varepsilon_t\)

Introduction

Quel estimateur utiliser pour estimer \(y_t = \beta_0 + \beta_1 x_t + \varepsilon_t\) ?

• Pour un PGD stationnaire, MCO ok.

• Si \(x_t\) et \(y_t\) sont DS, alors elles peuvent être cointégrées si c’est le cas –> MCE (modèle à correction d’erreurs, pour 2 séries uniquement). Si elles ne sont pas cointégrées (MCE impossible), mais quand même DS, alors on doit stationnariser les 2 séries et on emploie les MCO sur les séries stationnarisées.

Première étape : créer un chronogramme de la série étudiée

On cherche à détecter la présence (ou non) de 3 effets :

Tendance
Pour savoir si une série a une tendance, on compare le premier et le dernier point et on regarde si la tendance est haussière ou baissière. On regarde ensuite si la série fluctue autour d’une valeur, nulle ou non.

Effets saisonniers
On cherche à savoir si la série à des effets saisonniers (qui se reproduisent tous les ans). Cet effet doit se répéter plusieurs fois dans l’année et graphiquement doit se repérer par des pics réguliers. Si on en détecte, on les enlève pour la modélisation et on les remet à la fin.

Hétéroscédasticité
Est ce que la variance varie dans le temps ? On essaie d’insérer toutes les valeurs de la série entre 2 droites, et si ces dernières sont parallèles, alors la variance est finie (homoscedasticité).

On distingue 2 types de modèles

Modèle additif
\(x_t = m_t + s_t + u_t\)

Modèle multiplicatif
\(x_t = m_t s_t u_t\)

Avec :
- \(m_t\) la tendance
- \(s_t\) l’effet saisonnier
- \(u_t =\) un bruit où \(E(u_t) = 0\)

Exemple d’une série temporelle

Chronogramme de la température moyenne à Delhi entre 2013 et 2017

Décomposition de : l’effet saisonnier, la tendance et le bruit

Pour cela, on lisse notre série en calculant les moyennes mobiles puis en les retirant de cette dernière. Ensuite, l’effet saisonnier est calculée en faisant la moyenne, pour chaque unité de temps, de toutes les périodes. L”effet saisonnier est ensuite centrée. Enfin, le bruit est déterminée en supprimant la tendance et l’effet saisonnier de la série originale.

Le graphique de décomposition met en lumière la tendance haussière de la série sur la période donnée. On retrouve également un fort effet saisonnier sur les températures moyennes : des valeurs élevées en été et plus faibles en hiver. Les effets sont bien saisonniers étant donné que la fréquence des cycles est égale à la durée en année de la série : 4.

Simuler un AR, MA ou ARMA

bb = rnorm(150)
plot.ts(bb, col="darkblue")
abline(h=max(bb), col="red")
abline(h=mean(bb), col="green")
abline(h=min(bb), col="red")
abline(v=37, col="yellow")

acf(bb)

ts = arima.sim(list(order=c(1,0,0), ar=0.9), n=500)
acf(ts[200:500])

pacf(ts[200:500])

on test lunj box de l’ordre 1 à m, où m est le premier test significatif. On ne continue pas après car on aura tous les tests suivants significatifs = pas d’information supplémentaire

ts = rnorm(n=10000)
Box.test(ts, lag = 3, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  ts
## X-squared = 0.83205, df = 3, p-value = 0.8418

\[AR(2): f_t = 0.42f_{t-2} + \varepsilon_t\] \[ \varepsilon_t= - 0.42L^2 + 1 \] \[\Delta = b^2 - 4\times a \times c = 0 - 4 \times -0.42 = 1.68\] \[r_1 = \frac{-0 - \sqrt{1.68}}{2} \approx -0.65\] \[r_2 = \frac{-0 + \sqrt{1.68}}{2} \approx 0.65\] \[\mid z \mid = 0.65\]

df = ts(read.csv("serievendredi.csv"), freq=12, start=1950.1)
mean(df)

## [1] 0.1840205

sd(df)

## [1] 1.320298

plot.ts(df)

eacf(df)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o o o o o o o o  o  o  x 
## 1 x x x x x o o o o o o  o  o  x 
## 2 x x o o x o o o o o o  o  o  x 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x x x o o o o o o o  o  o  o 
## 5 o x x x x o o o o o o  o  o  o 
## 6 x x x o x o o o o o o  o  o  o 
## 7 x x x o x x x o o o o  o  o  o

x = acf(df)

x$acf[1:3]

## [1] 0.4166819 0.4573251 0.3223305

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

library(lmtest)

#ARMA avec p=2 et p=3, et m=1 et m=2 (fixed permet d'enlever l'ar1)
reg = Arima(df, order=c(3,0,3), fixed=c(0, NA, NA, NA, NA, NA, NA))
coeftest(reg)

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar2        0.330975   0.077722  4.2585 2.058e-05 ***
## ar3       -0.314124   0.077255 -4.0661 4.781e-05 ***
## ma1        0.117609   0.047284  2.4873   0.01287 *  
## ma2        0.294119   0.047706  6.1653 7.036e-10 ***
## ma3        0.819772   0.057920 14.1535 < 2.2e-16 ***
## intercept  0.173428   0.119980  1.4455   0.14832    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hist(df, breaks = 25, col = "lightblue")

plot(density(df))

Test sur les résidus

library(tseries)
jarque.bera.test(reg$residuals) #test de normalité

## 
##  Jarque Bera Test
## 
## data:  reg$residuals
## X-squared = 1.5301, df = 2, p-value = 0.4653

library(moments)

## 
## Attaching package: 'moments'

## The following objects are masked from 'package:TSA':
## 
##     kurtosis, skewness

agostino.test(reg$residuals)

## 
##  D'Agostino skewness test
## 
## data:  reg$residuals
## skew = -0.10868, z = -0.84582, p-value = 0.3977
## alternative hypothesis: data have a skewness

anscombe.test(reg$residuals)

## 
##  Anscombe-Glynn kurtosis test
## 
## data:  reg$residuals
## kurt = 3.2395, z = 1.0344, p-value = 0.301
## alternative hypothesis: kurtosis is not equal to 3

Box.test(reg$residuals, )

## 
##  Box-Pierce test
## 
## data:  reg$residuals
## X-squared = 0.55085, df = 1, p-value = 0.458

Si les résidus ne sont pas normalement distribués, les conclusions tirées sur la significativité des coef ne sont plus significatifs. IC à 95% autour des valeurs prévues. PLus ce dernier est petit, plus les prévisions sont précises.

SI autocorrélation (box test) à l’ordre 1, alors PAS BESOIN de faire des tests à un ordre supérieur (perte de points si incohérence). On doit alors refaire un test à un ordre plus élevé.

Test d’hétéroscédasticité conditionnelle de ARCH

Un arch(1) c’est quand la variance prévue à la date t dépend des aléas de t-1.

library(caschrono)
data(indbourse)
cac = as.numeric(indbourse[,5])
rt = diff(log(cac), rt=na.omit(rt))
plot.ts(rt)

library(FinTS)

## 
## Attaching package: 'FinTS'

## The following object is masked from 'package:forecast':
## 
##     Acf

ArchTest(rt, lags = 1)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  rt
## Chi-squared = 44.397, df = 1, p-value = 2.682e-11

A priori : pas d’effet arch sur les données annuelles car cela concerne plutot un effet qui sé déroule sur plusieurs jours. Meme chose que le box test : si effet arch à l’ordre 1, alors on en trouvera forcément à l’ordre 1+k.

Prévision

prev = forecast(reg, h=8, level=0.95)
plot(prev)

générer un processus TS + stationnarisrer

n = 1000
trend = seq(1,n)
xt = 3 + 1.2 * trend + rnorm(n=1000, sd=100)
plot.ts(xt[200:n])

reg = lm(xt~trend)
summary(reg)

## 
## Call:
## lm(formula = xt ~ trend)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -274.71  -70.70   -1.16   65.11  358.00 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.50240    6.39914  -0.704    0.482    
## trend        1.19936    0.01108 108.291   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 101.1 on 998 degrees of freedom
## Multiple R-squared:  0.9216, Adjusted R-squared:  0.9215 
## F-statistic: 1.173e+04 on 1 and 998 DF,  p-value: < 2.2e-16

xtstar = xt - reg$fitted.values
plot.ts(xtstar)

générer un processus DS (random walk) + stationnarisrer

RW <- function(N, x0, mu, variance) {
   z<-cumsum(rnorm(n=N, mean=0,sd=sqrt(variance)))
   t<-1:N
   x<-x0+t*mu+z
  return(x)
}
xt = RW(100, 0, 0.4,100)
plot.ts(xt)

Pour stationnariser un processus DS, on différencie à l’ordre 1 : DSstar = DSt - DSt-1. UN ds dépend de lui en -1 + un bruit.

dstart = diff(xt)
plot.ts(dstart)

Une série yt issu d’un pgd stationnaire est notée yt~I(0).

Une série xt issu d’un pgd DS est notée xt ~ I(1)

Normalement, après avoir stationnariser à l’ordre 1, on devrait avoir une série stationnaire.

Tests de racines unitaires

Méthode pour les racines unitaires :

1 - dickey-fuller lag=0, trend –> … –> voir arbre cours

Les tests de RU servent uniquement à savoir si le pgd générateur de nos données est ts, ds ou stationnaire.

\[X_t = \beta_0 + \beta_1 \times tendance_t + \rho X_{t-1} + \varepsilon_t\]

Pour que ce soit DS : \(\beta_1 = 0\) ET \(\rho=1\). R test si p-1 est différent de 0

TS si \(\beta_1\) significatif et \(\mid \rho \mid < 1\)

Pour savoir si beta1 est significatif, on regarde la p-value. Pour savoir si p-1 = 0, cette stat ne suit pas une student et il est donc interdit de regarder la p-value.

Dickey-fuller

library(caschrono)
library(urca)
data("indbourse")

cac=as.numeric(indbourse[,5])
rt=diff(log(cac))
rt=na.omit(rt)
summary(ur.df(rt, type = "trend", lags = 0))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.091507 -0.007683  0.000742  0.008310  0.099349 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.619e-05  1.013e-03   0.016    0.987    
## z.lag.1     -1.089e+00  2.996e-02 -36.361   <2e-16 ***
## tt          -8.839e-07  1.583e-06  -0.558    0.577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01685 on 1105 degrees of freedom
## Multiple R-squared:  0.5447, Adjusted R-squared:  0.5439 
## F-statistic: 661.1 on 2 and 1105 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -36.3614 440.7172 661.0757 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

beta1 pas significatif donc on test “drift”

summary(ur.df(rt, type = "drift", lags = 0))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.091637 -0.007579  0.000882  0.008334  0.099243 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0004738  0.0005062  -0.936    0.349    
## z.lag.1     -1.0891747  0.0299484 -36.368   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01684 on 1106 degrees of freedom
## Multiple R-squared:  0.5446, Adjusted R-squared:  0.5442 
## F-statistic:  1323 on 1 and 1106 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -36.3684 661.3313 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

beta0 pas significatif donc on test avec “none”

summary(ur.df(rt, type = "none", lags = 0))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.092131 -0.008043  0.000411  0.007860  0.098827 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## z.lag.1 -1.08846    0.02994  -36.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01684 on 1107 degrees of freedom
## Multiple R-squared:  0.5442, Adjusted R-squared:  0.5438 
## F-statistic:  1322 on 1 and 1107 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -36.3584 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Ici -36 < -1.95, donc pas RU dans les rt du CAC40.

On test sur une autre variable

cac=as.numeric(indbourse[,5])
rt=na.omit(cac)
summary(ur.df(rt, type = "trend", lags = 0))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -368.14  -33.98    2.61   37.02  366.62 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 39.426629  18.514694   2.129   0.0334 *
## z.lag.1     -0.006595   0.003114  -2.118   0.0344 *
## tt          -0.018665   0.009273  -2.013   0.0444 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67.18 on 1121 degrees of freedom
## Multiple R-squared:  0.004359,   Adjusted R-squared:  0.002582 
## F-statistic: 2.454 on 2 and 1121 DF,  p-value: 0.08643
## 
## 
## Value of test-statistic is: -2.1177 1.7381 2.4537 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

beta1 significatif. et présence de RU car -2.1>-3.4 (on accepte hypothèse nulle de présence de RU).

Cependant ces résultats ne sont valides que si les aléas de cette régression ne sont pas autocorrélés. Pour le déterminer, on fait acf et pacf.

plot(ur.df(rt, type = "trend", lags = 0))

Présence d’autocorrélation, donc on fait dickey-fuller augmenté: dicky-fuller avec des explicatives en plus (variable dépendante retardée jusqu’à plusieurs périodes).

\[X_t = \beta_0 + \beta_1 \times tendance_t + \rho X_{t-1} + \sum_{j=1}^p\gamma_j\Delta_{y_{t-j}} + \varepsilon_t\]

Première étape : déterminer p via MAIC (critère d’information valide). On commence par calculer q pour QDF pmax via la formule de Schwert.

data("indbourse")
cac=as.numeric(indbourse[,5])
rt=na.omit(cac)

pmax<-as.integer(12*(length(rt)/100)^(0.25))
library(CADFtest)

## Loading required package: dynlm

## Loading required package: sandwich

## Registered S3 methods overwritten by 'CADFtest':
##   method     from    
##   bread.mlm  sandwich
##   estfun.mlm sandwich

summary(CADFtest(rt,criterion="MAIC",type="drift",max.lag.y=pmax))

## Augmented DF test 
##                                              ADF test
## t-test statistic:                          -0.7020353
## p-value:                                    0.8441758
## Max lag of the diff. dependent variable:    5.0000000
## 
## Call:
## dynlm(formula = formula(model), start = obs.1, end = obs.T)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -374.37  -33.79    3.61   39.53  369.27 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.111601   9.702555   0.527  0.59842   
## L(y, 1)     -0.001463   0.002084  -0.702  0.84418   
## L(d(y), 1)  -0.098615   0.030183  -3.267  0.00112 **
## L(d(y), 2)  -0.038712   0.030295  -1.278  0.20158   
## L(d(y), 3)  -0.046830   0.030304  -1.545  0.12255   
## L(d(y), 4)   0.050961   0.030307   1.682  0.09295 . 
## L(d(y), 5)  -0.055033   0.030209  -1.822  0.06876 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67.18 on 1096 degrees of freedom
## Multiple R-squared:  0.02043,    Adjusted R-squared:  0.01507 
## F-statistic: 4.392 on 5 and 1096 DF,  p-value: 0.0005788

Ici il nous dit qu’il faut ajouter 5 variables

summary(ur.df(rt, type = "none", lags = 5))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -373.85  -33.59    4.51   38.65  369.82 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## z.lag.1     -0.0003661  0.0004295  -0.852  0.39420   
## z.diff.lag1 -0.0974478  0.0299202  -3.257  0.00116 **
## z.diff.lag2 -0.0386794  0.0300219  -1.288  0.19788   
## z.diff.lag3 -0.0469241  0.0300561  -1.561  0.11876   
## z.diff.lag4  0.0490590  0.0300758   1.631  0.10314   
## z.diff.lag5 -0.0556829  0.0299902  -1.857  0.06362 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66.87 on 1113 degrees of freedom
## Multiple R-squared:  0.01983,    Adjusted R-squared:  0.01454 
## F-statistic: 3.753 on 6 and 1113 DF,  p-value: 0.001055
## 
## 
## Value of test-statistic is: -0.8524 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

summary(ur.df(rt, type = "none", lags = pmax-1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -380.72  -34.29    3.25   38.32  379.35 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## z.lag.1      -0.0003698  0.0004352  -0.850  0.39565   
## z.diff.lag1  -0.0984377  0.0303759  -3.241  0.00123 **
## z.diff.lag2  -0.0336347  0.0304941  -1.103  0.27028   
## z.diff.lag3  -0.0453842  0.0305050  -1.488  0.13710   
## z.diff.lag4   0.0448165  0.0305307   1.468  0.14242   
## z.diff.lag5  -0.0559176  0.0305370  -1.831  0.06735 . 
## z.diff.lag6  -0.0159935  0.0305732  -0.523  0.60100   
## z.diff.lag7   0.0238343  0.0305712   0.780  0.43578   
## z.diff.lag8   0.0345928  0.0305979   1.131  0.25849   
## z.diff.lag9  -0.0473180  0.0306325  -1.545  0.12271   
## z.diff.lag10  0.0008217  0.0306656   0.027  0.97863   
## z.diff.lag11  0.0024622  0.0306680   0.080  0.93602   
## z.diff.lag12  0.0421615  0.0307345   1.372  0.17041   
## z.diff.lag13  0.0263415  0.0307444   0.857  0.39175   
## z.diff.lag14 -0.0027416  0.0307443  -0.089  0.92896   
## z.diff.lag15 -0.0256380  0.0307367  -0.834  0.40440   
## z.diff.lag16  0.0483451  0.0310409   1.557  0.11965   
## z.diff.lag17  0.0298639  0.0310894   0.961  0.33698   
## z.diff.lag18 -0.0297121  0.0310949  -0.956  0.33952   
## z.diff.lag19  0.0143978  0.0310938   0.463  0.64343   
## z.diff.lag20 -0.0198474  0.0310660  -0.639  0.52304   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67.19 on 1083 degrees of freedom
## Multiple R-squared:  0.03204,    Adjusted R-squared:  0.01327 
## F-statistic: 1.707 on 21 and 1083 DF,  p-value: 0.02438
## 
## 
## Value of test-statistic is: -0.8497 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

library(urca)
data("npext")
ts = ts(npext[,2], end=1988)
summary(ur.df(ts, lags = 0, type="drift"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.135530 -0.023892 -0.006165  0.018135  0.248059 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.066278   0.028331  -2.339  0.02089 * 
## z.lag.1      0.021740   0.007023   3.095  0.00242 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05508 on 126 degrees of freedom
## Multiple R-squared:  0.07067,    Adjusted R-squared:  0.0633 
## F-statistic: 9.582 on 1 and 126 DF,  p-value: 0.002422
## 
## 
## Value of test-statistic is: 3.0955 13.3256 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

On accepte l’hypothèse nulle de présence de RU.

pmax<-as.integer(12*(length(ts)/100)^(0.25))
library(CADFtest)
summary(CADFtest(ts,criterion="MAIC",type="drift",max.lag.y=pmax))

## Augmented DF test 
##                                             ADF test
## t-test statistic:                          0.8538915
## p-value:                                   0.9945694
## Max lag of the diff. dependent variable:   5.0000000
## 
## Call:
## dynlm(formula = formula(model), start = obs.1, end = obs.T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.208499 -0.010706 -0.000988  0.015055  0.092723 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.014502   0.024107  -0.602  0.54872    
## L(y, 1)      0.005377   0.006297   0.854  0.99457    
## L(d(y), 1)   0.777105   0.095496   8.138 7.11e-13 ***
## L(d(y), 2)  -0.355605   0.119722  -2.970  0.00366 ** 
## L(d(y), 3)   0.227033   0.122641   1.851  0.06685 .  
## L(d(y), 4)  -0.161553   0.119082  -1.357  0.17769    
## L(d(y), 5)   0.174038   0.096736   1.799  0.07477 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03666 on 109 degrees of freedom
## Multiple R-squared:  0.4856, Adjusted R-squared:  0.4573 
## F-statistic: 15.05 on 5 and 109 DF,  p-value: 3.215e-11

On refait alors dickey-fuller avec lag=2

summary(ur.df(ts, lags = 5, type="drift"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.206752 -0.012166  0.000132  0.014455  0.095946 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.022463   0.022830  -0.984  0.32718    
## z.lag.1      0.007267   0.005958   1.220  0.22505    
## z.diff.lag1  0.777168   0.086820   8.952 6.74e-15 ***
## z.diff.lag2 -0.316043   0.098307  -3.215  0.00169 ** 
## z.diff.lag3  0.207847   0.101924   2.039  0.04370 *  
## z.diff.lag4 -0.159571   0.097773  -1.632  0.10538    
## z.diff.lag5  0.115792   0.078547   1.474  0.14315    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03627 on 116 degrees of freedom
## Multiple R-squared:  0.5024, Adjusted R-squared:  0.4766 
## F-statistic: 19.52 on 6 and 116 DF,  p-value: 1.217e-15
## 
## 
## Value of test-statistic is: 1.2197 1.6976 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

On accepte la présence de RU (1.2197 > -2.88)

plot(ur.df(ts, lags = 5, type="drift"))

Comme les résidus sont autocorrélés, on diff.

ts = diff(ts)
summary(ur.df(ts, lags = 0, type="none"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.210762 -0.006813  0.004956  0.023481  0.139483 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## z.lag.1 -0.33336    0.06662  -5.004 1.84e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04526 on 126 degrees of freedom
## Multiple R-squared:  0.1658, Adjusted R-squared:  0.1592 
## F-statistic: 25.04 on 1 and 126 DF,  p-value: 1.844e-06
## 
## 
## Value of test-statistic is: -5.0042 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

On rejette l’hypothèse nulle de présence de racines unitaires –> stationnarité. On regarde si y’a de l’autoco dans les résidus

plot(ur.df(ts, lags = 0, type="none"))

Y’a toujours de l’autoco, on passe à ADF avec la spécification “none”.

pmax<-as.integer(12*(length(ts)/100)^(0.25))
library(CADFtest)
summary(CADFtest(ts,criterion="MAIC",type="none",max.lag.y=pmax))

## Augmented DF test 
##                                              ADF test
## t-test statistic:                          -2.4035396
## p-value:                                    0.0163226
## Max lag of the diff. dependent variable:    4.0000000
## 
## Call:
## dynlm(formula = formula(model), start = obs.1, end = obs.T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.218039 -0.004880  0.005037  0.021151  0.092171 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)   
## L(y, 1)    -0.20440    0.08504  -2.404  0.01632 * 
## L(d(y), 1)  0.02106    0.10632   0.198  0.84335   
## L(d(y), 2) -0.32332    0.10282  -3.144  0.00214 **
## L(d(y), 3) -0.06425    0.09448  -0.680  0.49793   
## L(d(y), 4) -0.21718    0.09303  -2.334  0.02139 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03699 on 110 degrees of freedom
## Multiple R-squared:  0.2422, Adjusted R-squared:  0.2078 
## F-statistic: 3.374 on 4 and 110 DF,  p-value: 0.01207

plot(ur.df(ts, lags = 4, type="none"))

summary(ur.df(ts, lags = 4, type="none"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.216926 -0.006945  0.005058  0.021283  0.093434 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## z.lag.1     -0.236927   0.080062  -2.959  0.00373 **
## z.diff.lag1  0.061736   0.087877   0.703  0.48373   
## z.diff.lag2 -0.243170   0.084092  -2.892  0.00456 **
## z.diff.lag3 -0.005375   0.076507  -0.070  0.94411   
## z.diff.lag4 -0.159021   0.074957  -2.121  0.03597 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03649 on 118 degrees of freedom
## Multiple R-squared:  0.2361, Adjusted R-squared:  0.2037 
## F-statistic: 7.295 on 5 and 118 DF,  p-value: 5.486e-06
## 
## 
## Value of test-statistic is: -2.9593 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Les économistes s’intéressent à la propriété de “mean reversion” (quand présent =TS = changement dans la tendance puis retour au régime initial) et à son absence (=DS = changement dans la tendance puis retour à une nouvelle tendance).

Spécification none : Ho=Ds, Ha=stationnaire (pas de mean reversion possible)

Spécification drift : Ho=DS, Ha=pas de tendance

Spécification trend : Ho=tendance linéaire même si \(\beta_1=0\), Ha=tendance linéaire même si \(\beta_1 \neq 0\) –> problème test permet une tendance linéaire sous Ha en incluant dans la régression la variable tendance alors que \(\beta_1\) pas pertinent sous Ho.

Dans DF, \(\beta_1\) et \(\beta_0\) sont des paramètres dits de nuisance. Dans DF, la stat \(t\) suit une loi normale et non-standard asymptotiquement sous Ho. Les distributions limites sont affectées par l’inclusions ou non de \(\beta_0\) et \(\beta_0 \times tendance\) (d’où des valeurs critiques différentes pour les 3 spécifications). Ces distributions sont fonction de processus de Wiener et les valeurs critiques sont obtenues par simulation.

L’intérêt de SP et LS est de s’affranchir de l’impact des paramètres de nuisance. Le test LS, le modèle à estimer (rho = 1 sous Ho) ne dépend pas de \(\beta_0\) et \(\beta_1\). Dans ce test, on a la possibilité de faire une estimation par bootstrap afin de nous affranchir de la distribution asymptotique. Bootstrap : ré-échantillonage avec remise (création d’échantillons à partir d’un échantillon initiale).

Exemple : si on a 35 individus avec un couple \(x_i - y_i\). Si on fait du MCO, on peut pas faire d’inférence car échantillon trop petit. Solution : on fait du bootstrap pour avoir une distribution des estimations.

Avantage de LS : on peut avoir 2 dates de rupture (par rapport à ZA).

Avec ZA : si on rejette Ho (pas de date de rupture sur la période), alors ça veut pas dire que y’a pas de RU. ça veut dire que en supposant que si y’a pas de date de rupture, alors y’a pas de RU. Dans la pratique : quasiment jamais de période sans changement structurel.

Faire la spécification break et crash, avec 1 et 2 dates pour les 2. On doit faire l’estimation par bootstrap. Entre crash et break, on prend le résultat que l’on a eu grâce à ZA (argument mymodel). L’argument mylags pas entre 1 et 5 (pmax est pas mal, mais il faudra probablement le diminuer car problème d’inversion de matrice)!. Le fait de chercher les dates de rupture peut être très dû au bootstrap.

#cac
library(caschrono)
data(indbourse)
cac = as.numeric(indbourse[,5])
cac = na.omit(cac)
rt = diff(log(cac))

#cpi
library(urca)
data("npext")
ts = ts(npext[,2])
rt = diff(ts)

#LS test
source("~/Desktop/M1S2/série temp/LeeStrazicichUnitRoot-master/LeeStrazicichUnitRootTest.R", encoding = 'UTF-8')
myBreaks <- 1
myModel <- "crash"
myLags <- 10
myLS_test <- ur.ls(y=rt, model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print")

## [1] -3.815446
## [1] "First possible structural break at position: 17"
## [1] "The location of the first break - lambda_1: 0.1 , with the number of total observations: 128"
## Critical values - Crash model:
##          1%     5%    10%
## [1,] -4.239 -3.566 -3.211
## [1] "Number of lags determined by general-to-specific lag selection: 10"
## Runtime:
## Time difference of 0.0008403818 mins

-3.8 < -3.5 : donc on rejette l’hypothèse nulle. si on fait print des résultats du test, on a les détails de l’estimation (avec lagmatrix = \(\rho\) et datmatNoLgas = \(\beta_1\))

Meme conclusions avec bootstrap

library(foreach)

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

library(doSNOW)

## Loading required package: iterators

## Loading required package: snow

library(parallel)

## 
## Attaching package: 'parallel'

## The following objects are masked from 'package:snow':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
##     parCapply, parLapply, parRapply, parSapply, splitIndices,
##     stopCluster

source("~/Desktop/M1S2/série temp/LeeStrazicichUnitRoot-master/LeeStrazicichUnitRootTestParallelization.R", encoding = 'UTF-8')
cl <- makeCluster(max(1, detectCores() - 1))
registerDoSNOW(cl)
myParallel_LS <- ur.ls.bootstrap(y=rt , model = myModel, breaks = myBreaks, lags = myLags, method = "Fixed",pn = 0.1, critval = "bootstrap", print.results = "print")

## [[1]]
## [1] -3.815446
## 
## [1] "First possible structural break at position: 17"
## [1] "The location of the first break - lambda_1: 0.1 , with the number of total observations: 128"
## Critical values - Crash model:
##          1%     5%    10%
## [1,] -4.239 -3.566 -3.211
## [1] "Number of lags used: 10"
## Runtime:
## Time difference of 0.002008883 mins

MCE

Si les 2 séries sont stationnaires, estimation par MCO valide.

Si les 2 sont TS, on transforme les 2 séries de telles sortes que les 2 séries soient stationnaires, puis on fait MCO sur les séries transformées (on retire la tendance \(y_t^* = y_t - reg\$fitted\)).

Si les 2 sont DS (les deux ont une tendance aléatoire), on se pose la question de est ce que les 2 séries partagent la même tendance ?

Si oui, alors on dit qu’elles sont cointégrées (elles partagent le même comportement de long terme, alors on doit estimer par MCE). Elles doivent être différenciées le même nombre de fois (généralement I(1) en économie) et il faut un vecteur cointégrant \(\beta\) stationnaire. On doit alors faire des tests de RU sur les résidus de l’estimation par MCO entre les variables (il ne doit pas y avoir de RU).

Si non, alors on ne fait pas d’estimation car fallacieuse.

Si les séries sont DS et TS, régression fallacieuse c’est mort.

Pour le projet 2, mettre les graphiques l’un au dessus de l’autre. Pensez à vérifier : cohérence dans les dates, les deux séries dans la même unité.